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FOREWORD 


This report has been prepared to expedite early domestic dissemination 
of the information generated under the contract. The data and conclusions 
must be considered preliminary and subject to change as further progress is 
made on this program. This is a progress report covering the work done 
during the first 12 months of the contract; it is not a final report. 

The NASA Program Manager is Dr. C.C. Chamis. 


ABSTRACT 


Accomplishments are described for the first year effort of a 5-year pro- 
gram to develop a methodology for coupled structural/therroal/electromagnetic 
analysis/tailoring of graded composite structures. These accomplishments 
include: (1) the results of the selective literature survey; (2) 8-, 16- , and 

20-noded isoparametric plate and shell elements; (3) large deformation struc- 
tural analysis; (4) Eigenanalysis ; (5) anisotropic heat tracsl^r analysis; and 
(6) anisotropic electromagnetic analysis. 
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1 . 0 INTRODUCTION 


This technical program is the work of the Engineering Mechanics and Life 
Management Section of the Aircraft Engine Business Group (AEBG) of the General 
Electric Company in response to NASA RFP 3-537260, "Coupled Structural/ 
Thermal/Electromagnetic (CSTEM) Analysis/Tailoring of Graded Composite Struc- 
tures." The overall objective of this program is to develop and verify analy- 
sis and tailoring capability for graded composite engine structures taking 
into account the coupling constraints imposed by mechanical, thermal, acous- 
tic, and electromagnetic loadings. 

The first problem that will be attacked is the development of plate and 
shell finite elements capable of accurately simulating the structural/thermal/ 
electromagnetic response of graded composite engine structures. Because of 
the wide diversity of engine structures and the magnitudes of the imposed 
loadings, the analysis of these is very difficult and demanding when they are 
composed of isotropic, homogeneous materials. The added complexity of direc- 
tional properties which can vary significantly through the thickness of the 
structure will challenge the state-of-the-art in finite element analysis. We 
are applying AEBG's 25 years of experience in developing and using structural 
analysis codes and the exceptional expertise of our University consultants 
toward the successful conclusion of this problem. To assist in. this, we are 
drawing heavily on previously funded NASA programs. 

We are drawing on NASA programs NAS3-23698, 3D Inelastic Analysis Methods 
For Hot Section Components, and NAS3-23687, Component Specific Modeling, in 
our development work on the plate and shell elements. In addition to these 
two programs, we will draw on NAS3-22767, Engine Structures Modeling Software 
System (ESMOSS) , and NAS3-23272, Burner Liner Thermal/Structural Load Model- 
ing, in Task III when we generate a total CSTEM Analysis System around these 
finite elements. This will guarantee that we are using the latest computer 
software technology and will produce an economical, flexible, easy to use 
system. 

In our development of a CSTEM tailoring system, we will build on NASA 
Program NAS3-22525, Structural Tailoring of Engine Blades (STAEBL) and AEBG 
program, Automatic Improvement of Design (AID) , in addition to the program 
system philosophy of ESMOSS. Because of the large number of significant 
parameters and design constraints, this tailoring system will be invaluable 
in promoting the use of graded composite structures. 

All during this program, we will avail ourselves of the experience and 
advice of our Low Observables Technology group. This will be particularly 
true in the Task V proof-of-concept . Their input will be used to assure the 
relevancy of the total program. 

Figure 1 shows our program and major contributions in flowchart form. 

This gives a visual presentation to the synergism that will exist between 
this program and other activities . 
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Figure 1. Program Flow Chart. 
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Figure 2 depicts an integrated analysis of composite structures currently 
under development in the composite users' community. The severe limitations 
of such a system are not highlighted because three major steps in the process 
are not shown. Figure 3 adds these steps. The analysis system really begins 
with a definition of geometry. A user then defines a finite element model 
simulating this geometry and the anticipated loading. The process then moves 
to defined Step 3. One cycle through the process ends with the prediction of 
individual ply average stresses and strains. Now comes a significant produc- 
tivity drain, namely, manual intervention to evaluate these stresses and 
strains against strength and durability limits. Based on this, the user must 
decide to (1) change the finite element model, (2) change the composite lami- 
nate, (3) both of the above, or (4) stop here. 

Obviously, there is a considerable cost savings to be obtained by select- 
ing Number 4. The CSTEM system will obviate the reasons for selecting Number 
4. This system, shown in Figure 4, begins with the definition of geometry, as 
before, but then proceeds to a definition of master regions which contain all 
of the necessary information about geometry, loading, and material properties. 
Step 3 is a constitutive model which develops the necessary structural, ther- 
mal, and electromagnetic properties based on a micromechanics approach. Fur- 
thermore, this constitutive model will contain the logic to generate the 
global finite element model based on the variation of the properties, as 
depicted in Figure 5. Using a nonlinear incremental technique, these global 
models will be solved for their structural, thermal, and ele^t. ©magnetic 
response. Based on this response the global characteristics will be evalu- 
ated, with convergence criteria and decisions made on remodeling. Once the 
global characteristics meet the accuracy requirements, the local characteris- 
tics are interrogated and decisions made on remodeling because of strength, 
durability, or hereditary effects. Once this cycle has been stabilized, opti- 
mization will be performed based on design constraint. Our goal in Task II is 
to develop finite elements whose characteristics make this system possible. 
Although the structural properties have been highlighted, the thermal and 
electromagnetic properties have as much or more variation, and less work has 
been done in these areas. 


1.1 EXECUTIVE SUMMARY 


Meetings were held with designers and management of the Low Observables 
Sections to determine their requirements. It was learned that the generic 
problem was that structures designed for optimum electromagnetic capabilities 
often have nonoptimum thermal and structural capabilities. There is a need 
for a tool that can accurately and efficiently analyze and iterate among 
these three fields so that viable compromise designs can be generated. At 
present, no such tool exists. 

Based on the Statement of Work and the results of the literature survey, 
we have established the 8-, 16- , and 20-noded isoparametric finite elements 
to be used to develop the CSTEM plate and shell element capabilities . These 
three elements meet the requirements for the plate elements. They have 
quadrilateral planform and are reducible to triangular planform, have nodes 
on the upper and lower surfaces, will meet accuracy requirements, and have 
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Figure 2. Composite Analysis System. 
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Figure 3. Total Composite Analysis System. 
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the proper degrees of freedom to perform all three types of analyses. The 
16- and 20-noded elements meet all of the requirements for the shell elements, 
including double curvature. 

Having established the basic building block, progress has been made in 
the three major areas, that is, structural, thermal, and electromagnetic, 
while maintaining the maximum commonality in computer software, such as shape 
functions, Jacobians, et cetera. 

Progress has occurred in the structural analysis area, as follows: 

• Defined the nonlinear equilibrium system of equations, including 
large deformation effects. Established an incremental, updated 
Lagranian solution. 

• Defined the overall programming architecture. 

• Wrote the modular stiffness matrix subroutines. 

• Wrote the modular mass matrix subroutines. 

• Wrote the subroutines for modular load. 

• Wrote subroutines for modular assembly and boundary conditions. 

• Established the input data format and wrote input subroutines. 

• Wrote the modular eigenvalue/eigenvector routines. 

• Incorporated the linear constraint equations. 

• Wrote a stress smoothing subroutine. 

• Wrote a modular equation solver. 

• Studied the micromechanics approach to stiffness formulation for 

composite elements. 

• Ran the verification and validation cases for the elasticity and 

eigenanalysis capabilities. 

• The file structure and data flow are currently under development. 

In the thermal area, the following progress has been made: 

• Anisotropic heat transfer equations have been defined for both 
linear and nonlinear conductivities and for steady-state and 
transient problems. 

• All of the subroutines necessary to perform a linear, steady-state, 
anisotropic heat transfer problem have been written, making maximum 
use of the code that is common with the structural analysis. 

In the electromagnetic area, the following progress has been made. 

• Results of the literature survey have been studied. 
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• Dr. M.V.K. Chari, General Electric 1 s expert on finite element elec- 
trical analysis, has been contacted. 

• Finite element equations for the electromagnetic field problem have 
been established. 
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2 . 0 TECHNICAL PROGRESS 


2* 1 TASK I - SELECTIVE LITERATURE SURVEY 

The first activity on this program was to perform a selective literature 
search using our internal General Electric data bases , the external COMPENDEX 
data base, and the literature supplied by the Texas A&M consultants. The 
pertinent articles turned up by this search are enumerated in Appendix A. 
Based on the results of this survey, it was proposed to the NASA Program 
Manager that the family of 8-, 16- , and 20-noded isoparametric elements be 
used for all three aspects of this program: structural, thermal, and electro 

magnetic. 

The proposed plan of attack was: 

• Work on three parallel efforts - structural, heat transfer, and 
electroma gnetic 

• Use an incremental updated Lagrangian approach for the large dis- 
placement structural problem. 

• Handle the coupling among the fields by an iterative procedure. 

2.2 TASK II - GRADED MATERIAL FINITE ELEMENTS 

2.2.1 Task IIA - Plate Elements 

The 8-, 16- , and 20-noded isoparametrics will be used as plate elements 
and shell elements. As plate elements, there will be a restriction that the 
midsurfaces be in a plane. This restriction will be the only difference 
between the plate and shell elements and, primarily, affects the program 
input. A simpler geometric and loading input can be used. Beyond this, all 
of the technical requirements are the same as for the shell elements. These 
will be discussed in the next section. 


2.2.2 Task IIB - Shell Elements 

AEBG has developed and used many different plate and shell elements. 

The elements to be used in this program are the 8-, 16- , and 20-noded iso- 
parametric elements. These have been used both as plates and as doubly curved 
shells for both linear and nonlinear material behavior. A review of these 
elements follows. 


Isoparametric Solid Elements 

The isoparametric solid elements permit the modeling of any general 
three-dimensional (3D) object, since the elements represent a discretization 
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of the object into finite elements which are 3D continuous representations. 

The basic term "isoparametric" means that the elements utilize the same inter- 
polating functions (also called "shape functions") to interpolate geometry, 
displacements, strains, and temperatures. It is, therefore, important that 
the user be aware that not just any displacement, geometry, and temperature 
field to be analyzed is necessarily compatible with a given element mesh. 

This is particularly true where high temperature or strain gradients occur. 

The following sections discuss the basic element formulation assumptions. 

Shape functions are used to describe the variation of some function G 
within an element in terms of the nodal point values. 

n 

G(x,y,z) = £ H i G i (x.,y i ,z i ) 
i=l 


where 

G(x,y,z) = the value of the function (such as displacement, tempera- 

ture) at any point with coordinates (x,y,z) within an 
element 

G (x ,y ,z ) = the value of the function at node point i 
i w i i 

H = the element "shape function" associated with node i 

i 

n = the number of nodes describing intraelement variation. 

In order to ensure monotonic convergence to the correct results, shape 
functions must satisfy several requirements. Satisfaction of these require- 
ments results in convergence from an upper bound. These displacement function 
requirements are: 

• They must include all possible rigid body displacements 

• They must be able to represent constant strain states 

• They must be differentiable within elements and compatible between 
adjacent elements. 

While the above conditions prove valuable for establishing upper bounds 
for solutions, they are not essential. Incompatible displacement modes are 
widely and successfully used. Their principal disadvantage is that stiffness 
may no longer be bounded from above. 

Curvilinear coordinates are introduced into the isoparametric concept to 
overcome the difficulty of formulating shape functions in global Cartesian 
coordinates. Also, generality in element geometry definition is obtained by 
this process. 
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A local curvilinear coordinate system (r,s,t), which ranges from -1 to 
1 within each element, is introduced in which shape functions are formulated. 
Also, a mapping from curvilinear to global coordinates is defined. A typical 
two dimensional element is shown in Figure 6. 

The same polynomial terms used in the Cartesian coordinates are used but 
with the curvilinear coordinates r,s,t replacing x, y, and z to generate 
shape functions . The r, s, and t coordinates are the same for all global 
element configurations . 


Typical finite element equilibrium equations: 

1. Structural 

[Mliiii) + [C]{u£> + [K]{ Ui } = {F B } + { F s } + {Fj} + { F c } + (F NL > 
[m] ■ / y p[H] T [H]dv "Consistent" mass matrix 

[C] ■ Damping matrix 

[K] " / v [B] t [D] [B] dv Stiffness matrix 

{F B > ■ / v t H l T t fB> dv Bod y forces 

{ Fg} ■ /s[H B ] T { fg}dS Surface tractions 
{ Ff} * / v [B] T [D]{ €T }dv Initial strains 

{FnL> “ / v 1 T {eNL) dv Nonlinear strains 
where 

{ ui}T = Iuj_ vj U 2 V 2 » 2 > — ] 

{u}T * [u v w] 

{u} * [H]{ui> 

{ e > ■ [B]{ui> 
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{a} - [D] {e e > - [D] ({e T0T }-{e THE ^ M }-{e N L}) 

< f B>T " U Bx *By fBzJ 

{ f S> T " t f Sx f Sy f SzJ 
Thermal 

[q]{Ti> + [K] {T£j - {Q B } + {Qg} + (Q C > 

[q] ■ /v c[H] T [H]dv Element heat capacity 

[K] - / v [B] t [R] [B] dv + / S2 o [H s ] T [H8]dS 2 

Element conductivity 

[Qfil " -f v rH] T {q B J dv Element heat generation 

(QsJ " / S 1 f HSl J T{ qs }dS l + f S «[H s ] T T B dS 2 

Boundary heat flow 

{Qc> * Concentrated heat flow 
where 

{ Ti }T - [ Tl t 2 t 3 ...] 

(T) * [HjiTj) 

(q) - [BjiTiJ 

C * Specific heat 

Tg * Boundary temperature 

oi ■ Surface heat transfer coefficient 
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Elect romagnetic 


YtSKAi) + t ID] (Ai) ♦ Y [E] (Aj) ♦ Jw [TKAi) - [T]{Ji> 

P 

where [S], ( D] , [E], and [T] are functions of [H]. 

{ > » Nodal point values of potential 

. Nodal point values of forcing functions 

Y * Reciprocal permeability 

w * Frequency 

p * Resistivity 

Other field problems have similar finite element expressions. 

8-Noded Solid 

The 8-noded solid element utilizes a formal mapping technique with dis- 
placement functions defined by a series of interpolating polynomials called 
"shape functions." In this manner, an arbitrary solid of six faces can be 
mapped from a parent coordinate system as shown in Figure 7. The numbering 
sequence must be as shown in this figure, except the user may define the nodes 
N1 and N5 as any convenient nodes on the solid. Note that node N1 has the 
parent coordinates (r,s,t) = (1,1,1). Once this has been established by the 
user, all face definition given later in this section is established. 

The shape functions for this element are as follows: 

H 2 = (1 + r) (1 + s) (1 + t)/8 

H 2 = (1 - r) (1 + s) (1 + t)/8 

H 3 = (1 - r) (1 - s) (1 + t)/8 

H 4 = (1 + r) (1 - s) (1 + t)/8 

H 5 = (1 + r) (1 + s) (1 - t)/8 

H 6 = (1 - r) (1 + s) (1 - t)/8 

H ? = (1 - r) (1 - s) (1 - t)/8 

Hg = ( 1 + r) (1 - s) (1 - t)/8 

Since this element has linear displacements on an edge, the ability to model 
the deformation due to bending is not possible, and the element will be overly 
stiff in bending. To overcome this, the addition of non-nodal degrees of 
freedom can be used as an option. These are called "incompatible modes," and 
are condensed out of the element stiffness after stiffness generation, leaving 
24 degrees of freedom per element. The interpolation for all element proper- 
ties and displacements is given by: ^ 
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“ = E Vi + (! - ' 2 ) a l + ( J - s2 ) a 2 + (' ’ t2 ) a 3 

i = 1 


8 

v = E H i V i + i 1 " r2 ) a 4 + { l ' s2 ) a 5 + i 1 ' ^ a 6 
i = 1 


= E H i w i + (l - f2 ) a 7 + ' g2 ) a 8 + 0 * t2 ) a 9 

i = 1 


where a through a are the extra "generalized” degrees of freedom which are 
condensed out. The user is given the option to include or not include the 
incompatible modes. In general, they should be included only when bending 
across an element is expected to be significant and only for elastic analyses. 

Given the coordinate system (r,s,t) as previously established, we can 
also now define the face numbering conventions and order of r^'es on a face. 
These definitions are needed to establish conventions for inputting pressure 
loads on the element and number of faces when displaying surface stresses on 
the faces. These conventions are summarized below: 


Face No. 

Location 

Nodes 

and Node 

Order 

on Face 

1 

r = +1 

N1 

N4 

N8 

N5 

2 

s = +1 

N1 

N5 

N6 

N2 

3 

t = +1 

N1 

N2 

N3 

N4 

4 

r = -1 

N7 

N3 

N2 

N6 

5 

s = -1 

N7 

N8 

N4 

N3 

6 

t = «1 

N7 

N6 

N5 

N8 

element has 

been formulated 

with 1 

variable 

temperature and general 


orthotropic material properties. During numerical integration for stiffness 
and equivalent nodal forces due to thermals, plasticity, and creep, the 
material properties at each integration point are evaluated at the temperature 
of that integration point. A Gauss integration scheme is used, and the user 
may choose an integration order of 2, 3, or 4 points in each direction 
(r , s , t) . 
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16-Noded Solid 


The 16-noded solid follows the same type of development for stiffness 
and loads as the 8-noded solid and, like the 8-noded solid, it has three dis- 
placement degrees of freedom per node, thus a total of 48 degrees of freedom 
on the element. Since this element has higher order interpolating functions 
in two directions on an edge, this element is sometimes called a "thick shell" 
element, as it can be used with reasonable approximation in this type of 
analysis where the shell thickness is in the t direction. 

The node numbering sequence must be as shown in Figure 8, except that 
the user can define the location of Nodes and N as desired. Note that 
Node N has the parent coordinates (r,s,t) - (1,1,1). Since this has been 
established, all face numbering is then defined. 

Starting with the basic interpolating functions for the corners, we 
define : 


G1 = (1 

+ r) 

(1 

+ s) 

(1 

+ 

t)/8 

G9 = (1 

+ r) 

(1 + 

s) 

(1 

- t)/8 

G3 = (1 

- r) 

(1 

+ s) 

(1 

+ 

t)/8 

Gil = (1 

- r) 

(1 + 

s) 

(1 

- t)/8 

G5 = (1 

- r) 

0 

- s) 

(1 

+ 

t)/8 

G13 = (1 

- r) 

(1 - 

s) 

Cl 

- t)/8 

G7 = (1 

+ r) 

(1 

- s) 

(1 

+ 

t)/8 

G15 = (1 

+ r) 

(1 * 

s) 

(1 

- t)/8 


Then, 

for the 

midside 

nodes 

i, the shape functions are: 

«2 = 

(1 - r) 

Cl 

+ s) 

Cl + 

t)/4 

«4 = 

(1 ~ r) 

(i 

- s) 

Cl + 

t)/4 

H 6 = 

(1 - r) 

Cl 

- s) 

C1 + 

t)/4 

»8 = 

(1 + r) 

Cl 

- s) 

C1 + 

t)/4 

t— ‘ 

o 

11 

(1 - r) 

C1 

+ s) 

C1 - 

t)/4 

»12 = 

(1 - r) 

C1 

- s) 

C1 - 

t)/4 

»14 = 

(1 - r) 

Cl 

- s) 

Cl - 

t)/4 

H 16 * 

(1 + r) 

Cl 

- s) 

o - 

t)/4 
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P84-73-4 



Figure 8. Sixteen-Noded Solid Coordinate and 
Node Numbering System. 
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and for the corner nodes, the modified shape functions are: 


H 1 =G 1 

- (H 2 * Hg)/2 


H 3 =G 3 

- (H 2 f H 4 )/2 


h 5 = g 5 

- <H 4 + H 6 )/2 


H 7 =G 7 

- (*6 + H 8 )/2 


«9 

- (H 10 ♦ H 16 )/2 


H 11 = G 11 

- (H 10 * H 12 )/2 


H 13 = G 13 

- (H, 2 ♦ H u )/2 

“I- 

H 15 " °15 

' < H 14 + "l6 )/2 

- = 


As in the case of the 8-noded solid, the user may use an option to 
specify that certain extra generalized degrees of freedom be added to intro- 
duce "incompatible modes" for bending. If through-thickness bending in the t 
direction is significant, these modes will prevent the element from being 
overly stiff. These generalized degrees of freedom are condensed out after 
stiffness formulation, leaving only the 48 nodal degrees of freedom. The 
complete element interpolating functions, used to interpolate displacements, 
etc . , are thus : 


u = 


16 

E 

i = 1 


H,U. + 
x x 


6 

E 

i = 1 


H . a . 
x x 


16 6 

v = Y H . V . + Y H.b, 

*S XX Is IX 

i = 1 i = 1 
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16 


6 


w = V H.W. + T H. c 

i—t X X 1 1 

i = 1 i = 1 

where the interpolating function coefficients for the generalized displace- 
merits are: 

Hj = r (1 - r 2 ) H 2 = s (1 - s 2 ) H 3 = (1 - t 2 ) 

H 4 = rs (1 - r 2 ) H 5 = rs (1 - s 2 ) H 6 = (1 - r 2 ) (1 - s 2 ) 

and the variables a 1 to a ^ , b 1 to b 6> and c ][ to are the 18 generalized * 
displacements whichare condensed out after stiffness formulation. 

Given the coordinate system (r,s,t) as previously established, we can 
also now define the face numbering conventions and order of nodes on a face. 
These definitions are needed to establish conventions for inputting pressure 
levels on the element and numbering of faces when displaying surface stresses 
on the faces. These conventions are summarized below: 


Face No. Location Nodes and Node Order on Face 


1 

r 

= 

+ 1 

N1 

N8 

N7 

N15 

N16 

hy 



2 

s 

= 

+ 1 

N1 

N9 

N10 

Nil 

N3 

N2 



3 

t 

= 

+ 1 

N13 

N2 

N3 

N4 

N5 

N6 

N7 

N8 

4 

r 

= 

-1 

N13 

N4 

N4 

N3 

Nil 

N12 



5 

s 

- 

-1 

N13 

N14 

N15 

N7 

N6 

N5 



6 

t 

= 

-1 

N13 

N12 

Nil 

N10 

N9 

N16 

N15 

N14 


This element has been formulated with variable temperature general ortho- 
tropic material properties. During numerical integration for stiffness and 
equivalent nodal forces due to thermals, plasticity, and creep, the material 
properties at each integration point are evaluated at the temperature of that 
integration point. A Gauss integration scheme is used, and the user may 
choose an integration order of 2, 3, or 4 points in each direction (r,s,t). 
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Tyenty-Noded Isoparametric Finite Element 


Figure 9 shows the representation and mapping of the 20-noded isopara- 
metric finite element in the the global (physical space) x, y, z coordinate 
system and the parent (barycentric) r, s, t, coordinate system. The trans- 
formation maps any global element into the same parent element with coordinate 
ranges of 

-1 < r < +1 

-1 < S < +1 

-1 < t < +1 

Figure 10 shows the node numbering in the parent coordinate system. Table 1 
gives the coordinate of each node. 


Let 


(u, v, w) = Displacements of a point in the (x, y, z) directions 


(u. , v., w. ) , i = 1 , 20 = Displacements of the nodes in the (x, y, z) 
111 directions 

{d} = [u, v, w] T 


(d. } = [u. , v. , w. ] 
i i’ l i 


Then for any point within the domain of the element: 


u 



= H. 


u. 


v 



H. v. * H. v. 

1 X XX 


w 



H. w . = H. w. 
i t 11 


u 


1 

o 

o 

•H 

» 

L 


1 

V 

s 

0 H. 0 

l 


v i 

w 

m m 


0 0 H. 

i J 


w, 

L ij 
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Table 1. Nodal Coordinates for 
20-Nodal Element. 


Node 

Number r s t 


1 

1 

-1 

-1 

2 

1 

1 

-1 

3 

1 

1 

1 

4 

1 

-1 

1 

5 

-1 

-1 

-1 

6 

-1 

1 

-1 

7 

-1 

1 

1 

8 

-1 

-1 

1 

9 

1 

0 

-1 

10 

1 

1 

0 

11 

1 

0 

1 

12 

1 

-1 

0 

13 

-1 

0 

-1 

14 

-1 

1 

0 

15 

-1 

0 

1 

16 

-1 

-1 

0 

17 

0 

-1 

-1 

18 

0 

1 

-1 

19 

0 

1 

1 

20 

0 

-1 

1 
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where the H^, in terms of the parent coordinate system are given as follows . 

The basic corner noded shape functions are 

Gj = (1+r) (l-s)(l-t)/8 
G 2 = (1+r) (1+s) (l-t)/8 
G 3 = (1+r) (1+s) (l+t)/8 
G 4 = (1+r) (1-s) ( 1+t ) / 8 
G 5 = (l-r)(l-s)(l-t)/8 
G 6 = (1-r) (1+s) (l-t)/8 
G ? = (1-r) (1+s) (1+t) /8 
G g = (1-r) (1-s) (l+t)/8 

The midside node shape functions are 
H 9 = ( 1+r) (l-s 2 )(l-t)/4 
H 1q = ( 1+r) ( 1+s) (l-t 2 )/4 
H X1 = (l+r)(l-s 2 )(l+t)/4 
H 12 = (1+r) (l-s)(l-t 2 )/4 
H 13 = (l-r)(l-s 2 )(l-t)/4 
H u = (1-r) (1+s) (l-t 2 )/4 
H 15 = (l-r)(l-s 2 ) (1+t)/ 4 
H lfi = (1-r) (1-s) (l-t 2 )/4 
H 1? = (l-r 2 )(l-s)(l-t)/4 
Hjg = (1-r 2 ) (1+s) (l-t)/4 
H 19 = (1-r 2 ) (1+s) (l+t)/4 
H 2q = (1-r 2 ) (1-s) (1+t )/4 


The modified corner node shape functions are 

H 1 = G 1 - <W H 17 )/2 
H 2 = G 2 - (H 9+ H 10+H)8 )/2 
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H 3 = G 3 - (H 10 « u « 19 )/2 

H 4 ~ °4 (“u +H 12 tH 20)/2 

H 5 * °5 - <W"l7 )/2 
«6 * C 6 - < H 13 +H 14* H J8 )/2 

"7 * °7 ' ' H 14 tH 15 tH 19 )/2 

K g = G g - (H 15 +H 16 +H 20 )/2 


In a three-dimensional context, the relations between the displacements 
and the strains are given by the following. 


Six Strain-Displacement Relations 
9u 

£ xx ” 9x 
_ 

E yy " 9y 

_ 9w 
£ zz ~ 3z 

_ i _ i / 3u + 3v\ 

£ xy 2 ^xy 2 y dy dx J 

1 _ 1 /8v 3w\ 

£ yz 2 ^yz 2 \3z 3y J 

1 1 /aw ^ 3u\ 

£ zx 2 ^zx ~ 2 \3x 3z / 


Therefore, the strain at any point within the elements domain is given 
by the following. 
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or 

{e} = [B.] {d.} 

e = B. d, 

1 1 

Since the Hi are defined in the parent coordinate system, we need to define 
the derivatives with respect to the global coordinate system by the chain 
rule of differentiation. 


3H. 

1 

3x 


3H , a 3H. * 3H . 

i 3r i 8s l 3t 

3r 3x 3s 3x 3t 3x 


3H. 

i 

3y 


3H . a 3H. a 3H. a . 

l 3r l 3s i 3t 

3r 3y 3s 3y 3t 3y 


3H . 

l 

3z 


3H. a 3H . a 3H, a . 

i 3r i 3s l 3t 

3r 3z 3s 3z 3t 3z 


To produce the required terms for the above expressions we need to develop 
the Jacobian matrix. 
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Inverting this Jacobian gives the following: 


9r 

3s 

at 

dx 

3x 

3x 

dr 

as 

3t 

dy 

3y 

3y 

dr 

3s 

at 

3z 

3z 

3z 


This is equivalent to the following: 


02 

at 


3z dz 
3s 3r 


J 


-1 


JL 

detJ 


3v 3y 3y 

9t 3s " 3r 


3x 3x 3x 

3t 3s 3r 


This gives all the necessary terms for the B matrix and therefore the ability 
to compute strain at any point within the element domain. 


e(x,y,z) = B(x,y,z). d ± 


e * 


3x 


0 0 


3H. 

i 

9y 


o o 


3H. 

i 

3z 


3H, 

2 

3x 


0 0 


0 0 


5 

ay 


0 0 


3H. 

i 

9z 


3H 


20 


2x 


0 0 


9H 

3y 


20 


0 0 


3H 

dz 


20 


20 


20 


J 20 
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The volume of the element is given by the following integral: 


V = 




detJ drdsdt 


The various structural, thermal, and electromagnetic quantities will be deter- 
mined by integrals over the volume, which takes the following form. 


F. . 
ij 





’ij 


detJ drdsdt 


Or, 


letting 


G. . 
ij 


g U 


detJ 


F. . 
ij 




G. . drdsdt 
ij 


This is normally evaluated by numerical integration of some form. 




EE 


G(r a’ 


V 


V 


W W u 
a b 


W 

c 


where m x n x o sampling points are used and W^, W c are the weighting 
factors for location r , s , t . One of our major efforts is the investi- 
gation of various numerical integration techniques with respect to the require- 
ments of CSTEM. Table 2 gives the locations and weights for up to the fourth 
order Glauss-Legendre Quadrature. Table 3 gives the information for up to 
sixth order Newton-Cotes (even interval) quadrature. 
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Table 2 • Gauss-Legendre Quadrature. 


Integration 


Order 

i 

r . 

w . 

— .J- 

2 

I 

I* 

+1A/3 

+1 


II 

-1/V3 

+1 

3 

I 

0 

8/9 


II 

+V0.6 

5/9 


III 

-Vo76 

5/9 

4 

I 

^3 + VO 

1 £o 

2 36 


II - 

y 3 + VO 

1 aM 

2 ~ 36 


III 

^ 3 - VA.i & 

1 + M 

2 36 


IV - 

^3 - VO 

i + M 

2 36 
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Table 3. Newton-Cotes Quadrature 


No. of 
Intervals 

n 


1 


2 


3 


4 


5 


6 


b 



n 



/ F< 

r)dr = 

(b-a 

>E 

C n F. 
1 1 

+ R 

n 

0 



i=0 



c“ 

c? 

c? 

Co 

c* 

C? 

0 

1 

2 

3 

4 

5 

1 

1 





2 

2 





1 

4 

1 




6 

6 

6 




1 

3 

3 

1 



8 

8 

8 

8 



7 

32 

12 

32 

_ 7 


90 

90 

90 

90 

90 


19 

75 

50 

50 

75 

19 

288 

288 

288 

288 

288 

288 

41 

216 

27 

272 

27 

216 

840 

840 

840 

840 

840 

840 


C 


n 

6 


41 

840 
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For the stiffness matrix: 

Kij »/ +1 f +l / +1 SiOBj detJ drdsdt 

where 

Gjj * B-j DBj detJ 

and D is the constitutive matrix 

The body force is given by: 

+1 +1 +1 T 

^B1 */ 2 / j / i^i F B detJ drdsdt 

where 

G i * h[ fB detJ 
and 

[fBl T * Cfflx F By F Bzl 

The initial strain effect is given by: 

F oi s / +1 / +1 / +1 B i° c o detJ drdsdt 
-1 -1 -1 

where 

r,j = Bi 0 c 0 detJ 


and 


[ £ 0 ] T - [Sc 'y h Y xy Y yz Y zxl 
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The "consistent" mass matrix is given by: 


m - r +1 r +1 T 

"iJ “ il -1 -1 pH H.detJdrdsdt 

* J 


Where 


G ij = pH i H j detJ 


The thermal strain effect is given by: 

+1 +1 +1 T 

F THi = h i l il BlDe TH detJdrdsot 


Where 


= B^De-^detJ 


And 


[ e THJ * [«« iT <x/ T c z iT 0 0 0 ] 


The nonlinear strain effect is given by: 


F NLi = {*’ fj' B T c NL detJdrsdt 

Where 


G i = B i De NL detJ 


And 


[€ NL ] " [ £ NLX e NLY e NLZ Y NLXY Y NLYZ Y NLZX ] 
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The surface traction effect is given by: 


sn 


■ f f "I. 


. f dA 
s 


-1 


Where the F . are the nodal forces due to surface tractions f on i 
sn s 

of constant r. Similar expressions apply for surfaces of constant e 


dA = 


Bx 


Bx 


By Bz. 

Bz By 

Bs 


at 


Bs Bt 

Bs Bt 

By 

X 

By 

dsdt = 

Bz Bx 

Bx Bz 

Bs 


Bt 

Bs Bt 

Bs Bt 

Bz 


Bz 


Bx By 

By Bx 

Bs 


Bt 


as at 

Bs Bt 

— - 


— - 






dsdt 


dA = Cdsdt 


Then 


Gi * H . f C 
Si s 


Large Deflection Theory 

Second order engineering strain-displacement equations 


Bu 1 
Bx 2 


3v 1 
By + 2 




3u\ 2 + ( §v\ 2 + / 3w^ 2 

3y/ \3y/ \9y . 


xy 


yz 


_ Bu + 3v + I 3u Bu 3v 3v Bw Bw I 
By Bx j Bx By Bx By Bx By | 

_ Bv | Bw 1 3u Bu | Bv Bv Bw 3w 

Bz By I < By Bz By Bz + By Bz I 

j 9u Bu Bv Bv Bw Bw [ 

Bx Bz Bx Bz Bx Bz 


Bw Bu 

e - — + — + 
zx Bx Bz 


surface 
or t . 
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Let 


[ , 1 ^ _ 9u 3v 3w 
x _0x 8x 9x _ 
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= [A] A { [G] [d] } = [A] [G] A [d] 

A [d 1 ] = [bJ A [d] 

A [e] = jj T jj A [d] 

M =[ B 0 + M 

|\] = [A] [G] 

The incremental stiffness is 

[R] = TkI + TkJ = f [B] T [D] [B] dV 


[ K J = JH m H dv 

V 

M ■ /(M T <°> M + M >»> M * f J T "» [ B ,]) d 

V X ' 1 

The force unbalance is given by: 
f(d) = j [B] [S] dV - [F] = 0 

V 

Where 

[F] are the applied forces 
A [el = [B] A [d] 
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A 



[S] dV + 


/ [B] T A [S] dV 
v 


But 


A [a] = [D] A [e] = [D] [B] A [d] 


And 


A [B] = A j^B 2 
Giving 


A [f] 


/ 4 M [SI dV + J [B] T [D] [B] a [d] dV 

V V 


Defining 


[k s ]a [di 


m [si dv 


Gives 


A [f] = ([ K s] * [ K l]) A Idl 

For eigenvalue problems: 

([ K i] * x [ K s]) a |d| * 0 

p 2 J T = [g) t [A] T 



A [d] = J [G] T 

V 


T 

A [A] 1 


[S] dv 


Then 



One of the major application areas for these graded material finite 
elements will involve small linear elastic strains but extremely large dis- 
placements. The severity of this condition is compounded by the involvement 
of temperature effects and body forces, conditions which have not been 
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extensively addressed. Additionally, eigenvalue/eigenvector analysis is 
required in the loaded conditions. We intend to take advantage of the sig- 
nificant improvements in computer speed and memory size in attacking this 
problem area. 

We are using an incremental-updated Lagrangian approach. The virtual 
work equation for the updated Lagrangian formulation leads to the nonlinear 
incremental displacements from time (t) to time (t + At). 


ijkl 5AE kl Ae ij ' 2 a ij 6 2Ae ik Ae kj 

- dV - 5 «< c +At > -/v t a u“ € ij d v 


where 


6w(t + At) =^ t < F 1 + AF i )5AU i dV 

+ / 8t (T 1 + AT 1 )5Au i ds 
l / SAUi 9AU 


Ae i j " 2 \ 3 x j ' 9 x i> 

f E Jki AAe kl^ £ i * s ec l u i va ^ ent to 

Jvt 1- ^ 3 

K] K! - (fn [ b l] t ^ M dV ) Kl 

-2 ! f ^ j 6Ae lk Ae kj dV is equivalent to 


(1) 


( 2 ) 


(3) 


(4) 
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CSd ( Auk l • (fit [\] T [°J [\] dv )l A " k } 

(5) 


f 0 , ,6AU, . -AU. ,dV is equivalent to 

#y£ IJ 



[ K N 2 ] M ■ (f vt [»2 ] T [° 2 ] L B 2 ] dV ) 1®"! 

(6) 


and 



f aj.<5Ae. .dV is equivalent to 

Jv t 1J J 


" 

f [\] T M dV 

(7) 

1 

J Vt 



where 



{4e( - [Bj {AU k l 

(8) 


Mil * [ b 2 ] f Auk| 

• 

(9) 


Mtj- [ AU l,X AU 1 ,2 AU 3,3] T 

(10) 


{au^} ** [ Au^ Au^ Au^...J T 

(11) 


The incremental equilibrium is 



([KJ + CS]) f Aok| ‘ f AP l 

(12) 


where 



is the linear stiffness matrix, i.e. 


- 

ia] -yjvT :<=] ca] dv 

(13) 


is the nonlinear stiffness matrix due to larger 

deformation, i.e. 

- 

IA] ‘ [Six] + IA 2 ] 

(14) 
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- 


( 15 ) 


PU ‘ / vt [\ f Pi] [\1 dV 
[ K N2J - / wt Plf [°2][ B 2] dv 

| Ap} is the incremental load from time (t) to (t + At) i.e. 

| Ap} = R(t + At) - J [Bj T |°j dV 


(16) 


The stress matrices are given by 


[Oo] 


[O] 


[ 0 ] 

0 

0 

[a] 

_ 0 

0 

- <*u 

a 12 

a 12 

a 22 

_ a i3 

a 23 


0 
0 

[O] 

a l 3 “| 

°23 


’33 J 


(17) 


(18) 


(19) 


foi] 


-2 a, 


11 


-2 a, 


22 


0 

0 

-2 a. 
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(Symmetric) 


~ °12 

0 


"°L3 

" a l2 

“°23 


0 

0 

" a 23 


"°13 

/ \ 

1 


1 

(°Ll +CT 22j 

2^3 


_ 2^3 


1 / \ 


1 

— 

2 ( a 22 + °33j 

1 

2 a 12 

/ , \ 



2 

( a 33 +0 llj 


( 20 ) 
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Eigenvalue/Eigenvector Analysis 


We are developing the algorithms for solving the lowest eigenvalue and 
eigenvector pairs for large structure system of equations; that is, 

K<|>. = A, M$. 

1 i x 

where 

K is the maximum symmetric stiffness matrix 

M is the maximum symmetric mass matrix 

is the ith eigenvalue 

$ is the ith eigenvector 

The methodologies coded so far are the determinant search and subspace itera- 
tion techniques . 

Both of these techniques have been investigated extensively in a finite 
element code that uses 8-noded shell elements and the COLSOL solution 
scheme. The addition of the eigenvalue/eigenvector analysis capability is 
closely linked to the solution procedure, and extensive changes to the COLSOL 
routines were required. Both lumped mass and consistent mass matrices were 
coded . 


Determinant Search Method 


In the determinant-search technique, a Sturm sequence check is very 
useful for estimating the proper eigenshift. Once the proper eigenshift is 
located in the neighborhood of the desired eigenvalue, a Rayleigh-quotient is 
then employed to iterate the trial vector until convergence is reached. A 
valid starting trial vector should be used, otherwise it may converge to the 
wrong result even though the proper eigenvalue is obtained in the repeated 
roots case. The logics for determining the proper eigenshift requires eigen 
deflation, Gram-Schmidt orthogonality if necessary, and schemes of accelera- 
tion, bisection, secant, or quadratic solutions. Appendix B shows the test 
case results of this procedure. 


Subspace Iteration 

The subspace iteration technique employs vector transformation into 
smaller q (n £ q 5 m) eigen spaces, that is 
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I * J 


L J 


The transformation matrix [x] is obtained through the successive inverse 
iteration process, 

[K] [x] k+1 = [m] [x] k 

where k is the iteration number. After transformation, the large eigenvalue 
problem is reduced to the small eigenvalue system containing the q eigenpairs 
and is easily solved by the well-known Jacobi rotation technique, that is 

[K] [v] = [m] [v] [A] 


where 

[K] = [x] T [K] [x] 

[m] = [x] T [m] [x] 

[x] = q x q diagonal eigenvalue matrix 

[v] = q x q eigenvector matrix corresponding [A] 

Even though these methodologies have been widely used, the solution accu- 
racy as well as economy rely heavily on the actual numerical implementation. 
The following numerical considerations are being kept in mind during program 
development . 


1 . Scale-Factors 

a. The mass value is, in general, much smaller than the stiffness 
value. A reasonable scale factor (10 6 for instance) should be 
imposed on the mass matrix in order to minimize numerical 
truncation errors. 

b. The determinant of the matrix [K - M] is, in general, very 
large. Care will be taken to avoid numerical overflow. 

2. Numerical Precision 


Depending on the computer hardware, double precision or even quad- 
ruple precision may be necessary. 
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3. 


Starting Trial Eigenvector 


The starting trial eigenvector significantly affects the iteration 
convergence rate. Most important is that the wrong trial eigenvec- 
tor may lead to the wrong solution in the determinant search tech- 
nique. In the subspace iteration technique, improper trial eigen- 
vectors may skip some eigenpairs. 

4. Convergence Criteria Limits 

a. Maximum number of iterations 

b. Rayleigh-Quotient iteration tolerance 

c. Eigenvector iteration tolerance 

d. Solution convergence error tolerance 

5. The eigenshift for determinant search Rayleigh-quotient iteration 
technique may converge to the wrong result if the proper eigenshift 
is not determined. In order to determine the proper eigenshift for 
the next eigenvalue, techniques should be investigated in order to 
minimize the number of large matrix decompositions. 

6. Efficient decomposition and back substitution matrix decomposition 
and back substitution play the most important role in the cost of 
performing eigenvalue analysis. 
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Heat Transfer Development 

Linear steady-state Equation: 
at time step: t + At: 



where 0 

K 

h 


90 9 © 96 1 
3x 3y 3z J 

K 0 0 

x 

0 K 0 

y 

0 0 K 

z 

convection coefficient 


s 


surface of the body 



governing equation: Linear heat transfer 


(K* 


K C ) 


<t+At 


= Q 


t+At 


<t+At 


(A) 


where 



the conductivity matrix 



K® B m dV 
m 


CD 


K 


c 


the convection matrix 



(II) 


the nodal point heat flow input vector 


Q (t+At) 


Q t+At Q g (t+At) + Q g (t+At) + Q c (t+At) 


(B) 


where Q = I 

gt+At m 


/ " 
J vm 


(m) . b(m) 

*t+At 


(III) 


Q g (t+At) 


I 

m 


/. 


S(rar S(m) 

(m) H a t+At 


ds 


(m) 


(IV) 


where Q “ a vector of concentrated nodal point heat flow input 

C 
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g 

The nodal point heat flow contribution Q t+ ^ t * s due to the convection 
boundary condition 


2 - the rate of heat generated in element 


s 

a 


- the surface heat flow input 


At that time: 


0 


( e ) _ 

t+At 


H (e) • 0 


t+At 


(la) 


© 


s(e) _ 
t+At 


H s( - e ' ) e 


t+At 


(lb) 


fl (e) _ R (e). f, 

e t+At “ 0 t+At 


(lc) 


where (e) — *. denotes element m 


0 


t+ ^ t — ► a vector of all nodal point temperatures at t+At 


0 


t+At 


0 ^t+At 0 2 t+ ^ t 0 ^t+At 


0 


m 


t+At 


] 


The matrix H 



Element Temperature 


B 



Temperature gradient interpolation matrice 


H 


s 



The surface temperature interpolation matrix 


51 


[B] — ► derivative of the shape function with respect to r, s, t and pre- 
multiplication by J 

- I f (m) “ ' 1 * 5< " )T • » S( " ) ' V tt4t > ’ ds( " ) (V) 

J sc 

© e t+ ^ t The given nodal point environmental temperatures, 
from this equation to find the Q^ +At (i.e. © e and h are given) 

Shape Function; 3 



H 1 = 8 1 ^ g 9 + 8 12 + *17^2 


H 2 ~ g 2 ^ g 9 + g 10 + * 18^2 


H 3 = g 3 ^ 8 10 + g ll + g 19^2 


H 4 = g 4 ^ g ll + g 12 + g 20^2 



H 5 = g 5 (g 13 + g 16 + g 17 )/ 2 
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H 6 = 8 6 ' (8 13 + 8 14 + 8 18 )/ 2 

H 7 = g 7 ^ g l4 + 8 15 + 8 19^2 

H 8 = g 8 ^ 8 15 + 8 l6 + 8 20^2 

H. = g. for ( i = 9.... 20) 

11 

g. = 0 if node i is not includes; otherwise: 
g i = G(r, r.^) G(s, s^ G(t, t^) 

GO, p^ = 1/2 (1+PiP) for P t = ± 1 

G(PO.) = (1-P 2 ) for p. = 0 


Jacobi’s operator [J] is 


The displacement functions: 


X = I H . x . 

l x 


Y = I H.y. 

i l 


Z = 1 H.z. 

i i 
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m 


§2 §2 

3* 3r 


9x 3^ §2 

38 38 


3x §2 

3t 3t 3t 



where r, s, t ate the local coordinates. 


[B] = [J] 1 


3M, 

3H. 

-a- ..... 

i 

3r 

3r 


Zl 

VU . 
1 

3s 

3s 

3H. 

3H. 

i 1 

i 

3t 

at 


the rate of heat generated per unit volume 

0 - the environmental temperature 

e 

£ S - the heat flow input to the surface of the body 
h - the convection coefficient 
q S = h(0 e - 0 S 3 


0 - is the temperature of the body 



— u( m ) Q 

®t+At " H ®t+At 


0 


s (m) _ H s(m) Q 


t+At 


0 -(m) = B (m). Q 


t+At 


0 , - a vector of all nodal point 

t+At 

temperatures at t+At 


0 


t+At 


= [Q : 


t+At 


05 


t+At 


a): from Equations (I), (II) we can find the 

(K k & K C ) 


(b): From Equation (B) the nodal point heat flow input 


^t+At ~ ^ fi (t+At) + ^g(t+At) + ^ C (t+At) 


where 


Q D (t+At) 


^ (t+At) 


= ljs2» 


vm H(m) *(t+At) dV 


s(m) s(m) .cOn) 

*(t+At) dS 


0 , A . - a vector of concentrated nodal point heat flow input 
c (t+At) 


It can be solved for 


(c): From Equation (V) the given nodal point environment temperatures 

®e(t+At) 
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^(t+At) " I 


sc 


h - H s(,n) • h 5 ^ 5 • e 


(t+At) 


dS 


(■) 


can solve fox 


e 

Q (t+At) 


(nodal point heat distribution) 


From items (a), (b) , (c), substitute into Equation A. 


(K 


+ K ) ®(t+At) = Q (t+At) + Q (t+At) 


The nodal point temperatures in each element can be found by using 
equations (la), (lb), (1c). 
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3D Magnetic-Field Problems - The general governing equations are: 


<J 

X 

II 

(1) 

i = pH 

(2) 

V • § = 0 

(3) 


where H, B, and 3 are vector quantities (with three space components); 
H = the magnetic-field strength, 

B = the magnetic-flux density, 


(Qnly 5 is a prescribed quantity; H and 
B represent six unknown variables.) 


, 

and J = the magnetic-current density; 
|j = the magnetic permeability; 


Define B = V x A, where t is a vector potential; then: 


V • B = V • (V x A) = 0 

VxH = Vx- = Vx-Vx^ = J 

M M 

V x i V x t = 3 (^) 

M 

This equation is equivalent to a stationarity of a variational functional II, 
defined as 


II 


-»T 1 

\ r A V X - v X a dv 

^ J v M 


- J A T 3 dv 

J V 


(5) 
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where t = .1- N,a. 

i=l i 1 

= shape function of x, y, z 

N t = N i (x, y, z) 

= the nodal value at Node i 

a, = a. (t) 
i i 

where v is all space. 


Define H = H + H where 
s m 


H = the magnetic source field, 
s 

an d = the induced magentizations . 


The H is any vector defined such that 
s 


V x H = 3 
s 

then Equation 1 becomes: 


( 6 ) 


V x H =0 
m 


(7) 


This can be identically satisfied by introducing a scalar potential <J> defined 
by: 


H = -V«D 


( 8 ) 


Eliminating B from Equations 2 and 3: 


V*§ = V'|j2 = V*|j(H + H ) = 0 
r r s m 

-V ♦ pV<|) + V • pS =0 


(9) 
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ORIGINAL PAQt IS 

OF POOR QUALITY 



For a simple scalar unknown (J> 


where 


Equation 


<t> = <l> (*• y, z. t) 

N i = N ± (x, y, z) 

*i = a i (t) 

9 may be written 


= §f (M x lb + af (M y iy } + af (M z H ) 


and, with the known quantity, written as: 


( 10 ) 


f = V • H • 3 (11) 

s 

For steady-state (time-dependent) problems, Equations 10 and 11 may be 
rewritten: 


§! <m x ff> + §y (M y I’ = £(x ' y ' z) 


( 12 ) 


z 



Boundary Conditions ~ Equation 12 is general and must be solved subject 
to additional constraints on the boundary surface. On S^, if nonlinear 
<t> = (x, y, z, t) , 
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4> = 4> (x, y, z) 


(13) 


while the remaining part of the boundary has the following condition on S 2 

H Q x + M y n y + H n z + 8 (x > z) = 0 (l4) 

where n , n , and n are the direction cosines of the outward normal of the 
surface, X and ^the union S x and S 2 forms the complete boundary T, (S x + S 2 = T). 
Equation 13 is called Dirichlet condition, and Equation 14 is called Cauchy 
boundary condition; if g = 0, it is called Neumann boundary condition. A 
field problem is said to have mixed boundary conditions when some portions of 
r have Dirichlet while others have Cauchy. 

In Equation 12, if (J = (J = (J = constant = p, V 2 <|> = f(x, y, z)/|j is 
called Poisson's equation. X If^ f = Z 0, then V 2 <|) = 0 is known as the Laplace 
equation. 

Variational Principal - The function <j>(x, y, z) that satisfies Equations 
12, 13, and 14 also minimizes the functional: 

1(40 = h J V [M X (|^) 2 + My(f^) 2 + *J Z (||) Z + 2f4>] dxdydz + /g 2 (g4)) dS 2 (15) 
6I(4») = 0 (16) 


Element Equations - Suppose that the solution domain v is divided into 
M elements of y nodes each. In each element, the unknown function may be 
expressed by 


4> (e) = Jj n. 4). = [N]{(|)} (e) 


07) 


where <{>. is the nodal value of (j) at Node i. Equation 17 implies that only 
nodal values of (J> are taken as nodal degrees of freedom, but derivatives of 
<}> may also be used as nodal parameters. 

n ( A 

(e) 

I(<}0 = I I (<}> ) 

e=l 


31 ( 4 > (e) ) = 0 

34). 

1 


i = 1, 2, y 
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For a Node i on boundary S 2 , from Equation 15 we have on Surface S 2 

91(<t> (e) ) _ 


(e) 


r r JL r^ (e \ + u ^ (e) r^ (e \ 

^(e) [M x 3x 3<)» i ( 3x 5 M y 3y a^/ay ^ P z 3z a<t> ± ( 3z 5 


+ f U. 1 dv(C) + 1 (e) (g U. 5 dSz 


(18) 


If Node i does not lie on S 2 , the second integral does not appear, 
to Equation 17, these terms typically become: 

a<fa (e) V 3N fe1 

I ■ ^ ST ♦! * t 8 "/ 8 *! <♦> 

1=1 


Referring 


_8 ,8* (e \ - 
3<|> ^3x J ~ 3x 


8t < * ) . N 

3<t>. "i 


Cel 

Thus on Surface (Equation 14) we have: 

ai(<i> (e) ) _ 

3«j>. 

1 


3N. 


3N. 


3N. 


J (e) [M x [ 3 N/ 3 x]{<|)}gji + My [3N/3y]{(|.}g^ + p z (3N/3z] + fN.] dv 


(e) 


+ J (e) (gNj dS^ 

S 2 


(e) 


Let <b = -I N.4». , then 
11 
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where 


[B(x , y, 2 )] = 


9Ni 

9N2 . - • • 



9x 

9x 

dx 

3Ni 

3N 2 

3N 

_K 

dy 

3y 

3y 

dNj 

3N 2 

3N 

_ 1 

3z 

3z 

dz 


From Equation 19, this may be rewritten 


where 


f [B] T [k ] [B] {<(>} dv + f [fN ] dv + J [gN ] dS 2 = 0 

V U Vi. 02 -L 



(nodal value) 


( 19 ) 


( 20 ) 
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0 


[k t ] = 


[k^.] = magnetic permeability on the principal axis (M x > Hy> an ^ M z ) • 


2.2.3 Task IIC - Graded Composite Materials 

An investigation into the merits of various numerical integration schemes 
has been pursued in order to determine which schemes might be best suited to 
the calculation of stiffness matrices for elements composed of several layers 
of different composite materials. All of the schemes investigated to this 
point integrate by summing weighted properties evaluated at sampling points 
which are spaced throughout the element volume. The differences between the 
various schemes is in the weights associated with the sampling points and the 
distribution of the sampling points. 

The modularization of the stiffness routines allows fairly easy implemen- 
tation of the various integration schemes. At present, the Gauss integration 
schemes and Newton-Cotes integration schemes have been coded. In addition, a 
selective Gauss integration scheme has been coded which uses different Gauss 
integration orders for calculation of normal stress/strain terms than the 
order used for shear stress/strain terms. Element stiffness matrices have 
been obtained using these methods, and a consistent, accurate method for com- 
paring them is being worked on. 

The previously coded stiffness routines for the 8-, 16- , and 20-noded 
isoparametric brick elements have been included in an existing finite-element 
code, and an extended checkout has begun. Elastic test case runs of isotropic 
materials have begun with comparisons to an existing finite-element code as 
well as the critical results. 


The test cases run to this point are from "A Proposed Standard Set of 
Problems to Test Finite Element Accuracy, " by R.H. MacNeal and R.L. Harder, a 
paper presented at the 25th SDM Finite Element Validation Forum, May 14, 1984, 
and are patch tests and cantilevered-element assemblages for eight-noded-brick 
elements (Figures 11 - 13), Tables 4 and 5 summarize the results of some of 
the test cases, and the theoretical results are presented in Table 6. It 
should be noted that all of these cases were run using single precision. 
Further improvement in the calculated answers is expected if double precision 
is used, and this will be investigated in the future. 

The 8-noded-brick results show that in the bending test cases the use of 
incompatible modes produces a considerably better solution as measured by tip 
displacement of the cantilevered beam model. However, the straightforward 
inclusion of incompatible modes causes the element to fail the patch test, as 
can be seen in the result for the existing code. An attempt to rectify this 
problem was made in the CSTEM code and is reflected in the results. Here, 
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•249 

.842 

.192 

.826 

•288 

.288 

.860 

.649 

.263 

.273 

.750 

.230 

.320 

.186 

.643 

.677 

.305 

.683 

.788 

.693 

.644 

.165 

.745 

.702 


Boundary Conditions: u • 10“ 3 (2* ♦ y ♦ i)/2 

* • 10~ 3 (i ♦ 2y ♦ t)/Z 
■ • 10~ 3 (« ♦ y ♦ 2z)/2 

Figure 11. Patch Test for Solids. 
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c) Parallelogrw Stupe Elements 


Length • €.0; Width ■ 0.2; Depth ■ 0.1 
E • 1.0 * 10 7 ; v - 0.30; %sh • 6 * 1 
Loading: Unit forces at free end 

Note: All elements have equal volume 


Figure 12. Straight Cantilever Beam. 





FIXED 




Inner Radius.* 4.12; Outer Radius • 4.32; Arc • 90* 
Thickness - 0.1; E * 1.0 x 10 7 ; v ■ 0.25; Nash *6x1 
Loading Unit forces at tip 


Figure 13. Curved Beam. 



Table 4. Existing FEM Code. 


Eight-Noded Bricks 


Without With 

Incompatibles Compatibles 


Patch Test 


0% 150% 

% Error in Stress 


Without With 

Incompatibles Compatibles 


Cantilevered Beam 



Extension 

0.9856 

(1.4%) 

0.9875 

(1.2%) 

Rectangular 

In Plane 

0.0928 

(90.7%) 

0.9922 

(0.8%) 

Elements 

Out of Plane 

0.0252 

(97.5%) 

1.070 

(7.0%) 


Twist 

0.8383 

(16.2%) 

0.7830 

(21.7%) 


Extension 

0.9845 

(1.5%) 

1.007 

(0.7%) 

Trapeze id 

In Plane 

0.0040 

(96.0%) 

0.1855 

(81.4%) 

Elements 

Out of Plane 

0.0067 

(99.3%) 

0.0286 

(97.1%) 


Twist 

0.6075 

(39.2%) 

0.6507 

(34.9%) 


Extension 

0.9846 

(1.5%) 

1.011 

(1.1%) 

Parallelogram 

In Plane 

0.0543 

(94.6%) 

0.7284 

(27.2%) 

Elements 

Out of Plane 

0.0084 

(99.2%) 

0.6367 

(36.3%) 


Twist 

0.4541 

(54.6%) 

0.7584 

(24.2%) 


Curved Beam 





In Plane 

0.0732 

(92.7%) 

0.9498 

(5.0%) 

Out of Plane 

0.2249 

(77.5%) 

0.7874 

(21.3%) 


Normalized Tip Displacement (% Error) 
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Table 5. CSTEM Stiffness Code. 


Eight-Noded Bricks 



Without 

With 


Incompatibles 

Compatibles 

Patch test 

0% 

0% 


% Error 

in Stress 




Without 

With 



Incompatibles 

Compatibles 

Cantilevered Beam 




Extension 

0.9856 (1.4%) 

0.9875 (1.2%) 

Rectangular 

In Plane 

0.0928 (90.7%) 

0.9922 (1.8%) 

Elements 

Out of Plane 

0.0252 (97.5%) 

0.9081 (9.2%) 


Twist 

0.8394 (16.1%) 

0.8186 (18.1%) 


Extension 

0.9845 (1.5%) 

1.007 (0.7%) 

Trapezoid 

In Plane 

0.0040 (96.0%) 

0.1851 (81.5%) 

Elements 

Out of Plane 

0.0067 (99.3%) 

0.0285 (97.1%) 


Twist 

0.6078 (39.2%) 

0.6502 (35.0%) 


Extension 

0.9846 (1.5%) 

1.009 (0.9%) 

Parallelogram 

In Plane 

0.0543 (94.6%) 

0.7343 (26.6%) 

Elements 

Out of Plane 

0.0084 (99.2%) 

0.5960 (40.4%) 


Twist 

0.4543 (54.6%) 

0.8008 (19.9%) 


Curved Beam 





In Plane 

0.0733 

(92.6%) 

0.9520 

(4.8%) 

Out of Plane 

0 . 2265 

(77.4%) 

0.8771 

(12.3%) 


Normalized Tip Displacement (% Error) 
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Table 6. Theoretical Solutions. 


Patch Test 

ex = ey 
ax = cry 


Straight Beam 

Extension 
In Plane 
Out of Plane 
Twist 


Curved Beam 

In Plane 
Out of Plane 


ez = yxy = yyz = yzx = 10 3 

az = 2000 psi, txy = xyz = Tzx = 400 ps 


(Tip Displacement in Direction of Load) 

3.0 x 10 -5 in. 

0.1081 in. 

0.4321 in. 

0.00341 in. 


(Tip Displacement in Direction of Load) 

0.08734 in. 

0.5022 in. 


the incompatible modes are included in the stiffness matrix and thus in the 
solution of the system of equations. This improves displacements; however, 
the recovery of stress from the displacement solution doesn T t include 
incompatible modes, so the patch test is satisfactory. Along with other 
possibilities, this technique is still being investigated. 

Another observation that can be made is the reinforcement of the fact 
that modeling technique can greatly influence the results. This is evident 
from the use of trapezoid- and parallelogram-shaped elements to solve the 
same cantilever beam problem. A trapezoid-shaped element in particular 
should be avoided in view of the bad results. 

The 20-noded brick passes the patch test with no problem regardless of 
integration order or whether single or double precision is used. The 20-noded 
brick is a quadratic displacement element by nature, and so there are no 
incompatible modes associated with it. The problem with the patch test as 
exhibited in the 8-noded brick is then avoided by the 20-noded brick. 

The regular beam was used for the remainder of the investigation. 

Because of the single layer of elements, reduced integration produces poor 
results when looking at displacements. This is due to the presence of zero 
energy modes: spurious nodal displacements which still satisfy the elemental 

equations at the Gauss points. More than one layer of elements would help to 
eliminate this phenomena due to the continuity of the element boundaries. 

It was found that the use of single or double precision would greatly 
affect the results as can be seen in the tabulated results, Table 7. This 
points to a numerical sensitivity to rounding off or truncation. Using the 
CSTEM code, attempts were made to improve the results by performing certain 
operations in double precision while the majority of operations remained in 
single precision. It was finally found that the best improvement was obtained 
by calculating the shape functions using double precision while all other 
operations remained single precision. 


2.2.4 Task IIP - I/O and Solution Techniques 

To achieve computational efficiency and to avoid repeated codings, a 
modular equation solver has been written to perform combinations of the fol- 
lowing functions: 

• Matrix decomposition (matrix triangularization) 

• Force vector reduction and back-substitution 

• Matrix decomposition and force vector back-substitution 

• Partial matrix static condensation 

• Matrix determinant and Sturm sequence count. 

Investigations of the eigenvalue/eigenvector large deformation structural 
problem were conducted. Two popular techniques were studied: determinant 

search and subspace iteration. 
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EXISTING FEM CODE 
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Out of Plane 0.3030(69.7%) 0.6043(40.0%) 0.9609(3.9%) 0.9352(6.5%) 

Twist 0.5742(42.6%) 0.8406(15.9%) 0.8516(14.8%) 0.5363(46.4%) 



Determinant Search - This technique employs inverse iteration in conjunc- 
tion with eigenshift and Rayleigh quotient method. This method relies heavily 
on the equations solver to compute Sturm sequence checks and matrix determi- 
nants. The starting trial eigenvector plays a very important role in converg- 
ing to the appropriate corresponding eigenpair. Orthogonality criteria must 
be imposed for computing repeated, or cluster, eigenpairs. 

Subspace Iteration - This method employs inverse iteration, subspace 
transformation, and Jacoby successive rotation. Besides the efficient equa- 
tions solver, this method relies heavily on the stiff/mass ratio and conver- 
gence criteria. This method will generally compute the desired eigenpairs, 
but it may not guarantee the lowest eigenpairs. 

Constraint equations play a very important role in solving boundary value 
problems. Two techniques have been developed for solving FEM structural prob- 
lems with constraints. The first technique is the so-called "penalty func- 
tion" method which treats the constraint equation as a stiff finite element 
and is very easy to incorporate in the program, such as: 

+ To construct a fairly stiff symmetric stiffness matrix out of the 
given linear constraint equation 

• To assemble the constraint stiffness matrix into the global overall 
structural stiffness matrix. 

The technique generally yields satisfactory solutions as long as the 
number of constraints is not large. 

The second technique is the classical matrix partitioning method. For 
computational efficiency, Gauss elimination is used to perform static con- 
densation rather than using matrix inversion to effect matrix partitioning 
and reduction. This method is mathematically exact and should be employed 
when the number of constraints becomes large. However, implementation of 
this technique into the program may be complicated due to the variation of 
structural modelings. 

For those stress computations within the finite element, Lagrange inter- 
polation demonstrated very good results when data are interpolated or extrapo- 
lated from the known Gaussian quadrature quantities. 

In order to extend the program capability from static analysis to eigen- 
value/eigenvector dynamic analysis, the assembly of consistent mass matrices 
is required. Because of the similarity between element stiffness matrix and 
element consistent mass matrix * the assembly routines were extended to accom- 
modate either. 

A new data file structure has also been added to the CSTEM stiffness 
routines. The investigation into different integration techniques emphasized 
the need to base storage of certain data on the integration point, rather than 
the element as previously done. This is due to the fact that the number of 
integration points may vary from element to element, depending on the inte- 
gration scheme used. The new data-file structure was incorporated into the 
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code used for the preceding investigations, so installation can be considered 
complete . 


2.2.5 Task HE - Stand-Alone Codes 

These capabilities are being generated as stand-alone codes while the 
development continues. The major effort in this area is the development of 
the file structure and data flow for this complex, interconnected problem. A 
change in the file structure is under consideration as a result of the larger 
number of integration points being considered. The present element-based 
file structure will not adequately handle this problem. 
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APPENDIX B - EIGENVALUE SOLUTIONS 


Closed-form eigenvalue solutions of cantilever continuous beam with 
rectangular cross section 



1 


w 2 

1 



El 

Y 


where 

E = Young f s modulus 

3 

I = Area moment of inertia = wt / 12 
\ - Mass per unit length = pwt 

p = Mass density 

a. = Constants = 1.875, 4.694, 7.855, ■ 


1 . 

2 . 


Finite Element Models and Results 

Four MSS8 shell elements of equal length were used. 
Material data: 

E = 30 x 10 6 psi 


p = 0.298/386 


lb-sec' 

, 4 
in 
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3. 


Constant cantilever dimensions: 


L = 8 in. 

W = 2 in. 

4. Results: 

Test Case Thickness 

t = 2 in. 
t = 1 in. 

t = 0.2 in. 


Theoretical 
= 39 

A = 1535 

Aj = 9.77 
\ 2 = 384 

= 0.39 
A 2 = 15.4 

k 3 = 120 
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MSS8 

\ 1 = 38.9 
\ 2 = 1014 

Aj = 10.3 
A 2 = 347 

Aj = 0.42 
A 2 = 15.7 
X 3 = 118 
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